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Abstract 

We address the thermodynamics (equilibrium density profiles, phase diagram, insta- 
bility analysis...) and the collapse of a self-gravitating gas of Brownian particles in D 
dimensions, in both canonical and microcanonical ensembles. In the canonical ensemble, 
we derive the analytic form of the density scaling profile which decays as f(x) ~ x~ a , 
with a = 2. In the microcanonical ensemble, we show that / decays as f(x) ~ x _Qmax , 
where a max is a non-trivial exponent. We derive exact expansions for o max and / in the 
limit of large D. Finally, we solve the problem in D = 2, which displays rather rich and 
peculiar features. 



1 Introduction 

In an earlier paper we studied a model of self-gravitating Brownian particles confined within 
a three-dimensional spherical box. We considered a high friction limit in which the equations 
of the problem reduce to a Smoluchowski-Poisson system with appropriate constraints ensuring 
the conservation of energy (in the microcanonical ensemble) or temperature (in the canonical 
ensemble) B. The equilibrium states (maximum entropy states) correspond to isothermal 
configurations which are known to exist only above a critical energy or a critical temperature 
(see, e.g., ||). When no hydrostatic equilibrium exists, we found that the system generates 
a finite time singularity (i.e., the central density becomes infinite in a finite time) and we 
derived self-similar solutions describing the collapse. This study was performed both in the 
microcanonical and canonical ensembles, with emphasize on the inequivalence of ensembles for 
such a nonextensive system. In the canonical ensemble, we showed that the scaling exponent 
for the density is a = 2 and we determined the invariant profile f{x), satisfying f(x) ~ x~ a for 
x — > +oo, analytically. In the microcanonical ensemble, the scaling exponent a ~ 2.21... and 
the corresponding invariant profile f(x) were determined numerically. 

In this paper, we propose to extend our analysis to a space of arbitrary dimension D. 
The interest of this extension is twofold. First, we shall consider an infinite dimension limit 
D — > +oo in which the problem can be solved analytically. In particular, it is possible to 
determine the scaling exponent a(D) and the profile f(x,D) in the microcanonical ensemble 
by a systematic expansion procedure in powers of D~ x (Sec. |3.4j ), while the canonical value is 
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always a = 2 and the profile can be calculated exactly for any dimension (Sec. |3.3| ). We show 
that, already up to order 0(D~ 2 ), the results of the large D expansion agree remarkably well 
with those found numerically for D = 3. Moreover, we show that the nature of the problem 
changes at two particular dimensions D = 2 and D = 10. In Sec. |^, we compute the equilibrium 
phase diagram as a function of the dimension. For 2 < D < 10, the T — E curve has a spiral 
shape like in 3D. For D > 10 and D < 2, the T — E curve is monotonous. The dimension 
D = 2 is critical and requires a particular attention that is given in Sec. f|. We show that 
for D = 2 the system generates a Dirac peak (containing a finite fraction of mass) while for 
D > 2, the central singularity contains no mass (at the collapse time). The case D = 2 has 
interest in theoretical physics regarding 2D gravity [|J and string theory || (in connection with 
the Liouville field theory). It has also applications in the physics of random surfaces || and 
random potentials 0, 2D turbulence || and chemotaxis (for bacterial populations). Finally, 
the dynamical equations considered in this paper and in Ref. |§ are receiving a growing interest 
from mathematicians who established rigorous results concerning the existence and unicity of 
solutions for arbitrary domain shape without specific symmetry. We refer to the papers of 
Rosier []IIJ and Biler et al. [II]], and references therein, for the connection of our study with 



mathematical results. 



2 Equilibrium structure of isothermal spheres in dimen- 
sion D 

2.1 The maximum entropy principle 

Consider a system of particles with mass m interacting via Newtonian gravity in a space of 
dimension D. The particles are enclosed within a box of radius R so as to prevent evaporation 
and make a statistical approach rigorous. Let /(r, v, t) denote the distribution function of the 
system, i.e. f(r,v,t)d D rd D v gives the mass of particles whose position and velocity are in the 
cell (r, v; r + d D r, v + d D v) at time t. The integral of / over the velocity determines the spatial 
density 



(1) P = Jfd D v. 
The total mass of the configuration is 

(2) M = f pd D r. 



In the meanfield approximation, the total energy of the system can be expressed as 

(3) E = \j dDrd ° W + 2 / P$ = K + W ' 

where K is the kinetic energy and W the potential energy. The gravitational potential $ is 
related to the density by the Newton-Poisson equation 

(4) A$ = S D Gp, 

where is the surface of a unit sphere in a D-dimensional space and G is the constant of 
gravity. Finally, we introduce the Boltzmann entropy 

(5) S = - [ f In fd D vd D v, 
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and the free energy (more precisely the Massieu function) 

(6) J = S — pE, 

where (3 = l/T is the inverse temperature. If the system is isolated, the equilibrium state 
maximizes the entropy S at fixed energy E and mass M (microcanonical description). Alterna- 
tively, if the system is in contact with a heat bath which maintains its temperature fixed, the 
equilibrium state maximizes the free energy J at fixed mass M and temperature T (canonical 
description). It can be shown that for systems interacting via a long-range potential like gravity, 
this meanfield description is exact so that our "thermodynamical approach" is rigorous. 

To solve this variational problem, we shall proceed in two steps. We first maximize S (resp. 
J) at fixed M, E (resp. T) and p(r). This yields the Maxwell distribution 

(7) /= (2^75 p(r)c "^- 

It is now possible to express the energy and the entropy in terms of p(r) as 

(8) E = ^MT + 1 



J p<t>d D r, 



(9) S = ^-MlnT - J p\npd D r, 

where we have omitted unimportant constant terms in the entropy @. The entropy and the 
free energy are now functionals of p(r) and we consider their maximization at fixed energy or 
temperature. Introducing Lagrange multipliers to satisfy the constraints, the critical points of 
S or J are given by the Boltzmann distribution 

(10) p = Ae~^. 

Then, the equilibrium state is obtained by solving the Boltzmann-Poisson equation 

(11) A$ = S D GAe~^, 

and relating the Lagrange multipliers to the appropriate constraints. Note that a similar 
variational problem occurs in the context of two-dimensional turbulence (D = 2) to characterize 



large-scale vortices considered as maximum entropy structures |TJ, ||, [U 

It is easy to show that there is no global maximum of entropy at fixed mass and energy 
for D > 2 (see, e.g., Ref. [0J for D — 3). We can make the entropy diverge to +oo by 
approaching an arbitrarily small fraction of particles in the core (M core <C M) so that the 
potential energy goes to — oo. Since the total energy is conserved, the temperature must rise to 
+oo and this leads to a divergence of the entropy to +oo. Note that if we collapse all particles 
in the core, the entropy would diverge to — oo. Therefore, the formation of a Dirac peak is 
not thermodynamically favorable in the microcanonical ensemble. For D = 2, there exists a 
global entropy maximum for all energies. On the other hand, there is no global maximum of 



free energy at fixed mass and temperature for D > 2 (see, e.g., Ref. |18j for D = 3) and if 
T < T c = GM/4 for D = 2 (see Appendix 0). We can make the free energy J diverge to 
+oo by collapsing all particles at r = 0. Therefore, a canonical system is expected to form a 
Dirac peak. For D = 2 and T > T c , there exists a global maximum of free energy. For D < 2, 
there exists a global maximum of entropy and free energy for all accessible values of energy and 
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temperature. We refer to Refs. [14], [15] for a rigorous proof of these results. When no global 



maxima of entropy or free energy exist, we can nevertheless look for local maxima since they 
correspond to metastable states which can be relevant for the considered time scales. Of course, 
the critical points of entropy at fixed E and M are the same as the critical point of free energy 
at fixed T and M. Only the onset of instability (regarding the second order variations of 5* 
or J with appropriate constraints) will differ from one ensemble to the other. For D = 3, this 
stability problem was considered by Antonov |16[ and Padmanabhan |T7| in the microcanonical 
ensemble and by Chavanis |18| in the canonical ensemble, by solving an eigenvalue equation. It 
was also studied by Lynden-Bell & Wood |H| and Katz |20] by using an extension of Poincare 
theory of linear series of equilibria. We shall give the generalization of these results in Sec. [2.6| 
to the case of a system of arbitrary dimension D. 



2.2 The D-dimensional Emden equation 



To determine the structure of isothermal spheres, we introduce the function ip = /3($ — $o)j 
where <£>o is the gravitational potential at r = 0. Then, the density field can be written 



12) 



P = Pod 



where po is the central density. Introducing the notation £ = {SDGf3po) l / 2 r and restricting 
ourselves to spherically symmetrical configurations (which maximize the entropy for a non- 
rotating system), the Boltzmann-Poisson Eq. (|TTD reduces to the form 



(13) 



d 



e- 1 # 

which is the Z)-dimensional generalization of the Emden equation |2T 
has a simple explicit solution, the singular sphere 

. _ 2(D - 2) 



For D > 2, Eq. (TTJ) 



(14) 



e 



The regular solution of Eq. ( |1~3"D satisfying the boundary conditions 
(15) V = = at £ = 0, 

0, we can expand the solution in Taylor series and we 



must be computed numerically. For £ 
find that 



(16) 



if, 
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2ZT 8D(D + 2 



D + l 



To obtain the asymptotic behavior of the solutions for £ 
t = ln£, ip = 2 ln£ — z changes Eq. ([13]) in 



24 D 2 (D + 2)(D + A) 

+oo, we note that the transformation 



(17) 



^ + (Z>- 2)^ + 2(1, -2). 



For D > 2, this corresponds to the damped oscillations of a fictitious particle in a potential 
V(z) = e z — 2{D — 2)z, where z plays the role of the position and t the role of time. For 
t — > +oo, the particle will come at rest at the bottom of the well at position z = ln(2(D — 2)). 
Returning to original variables, we find that 



(18) 



2(D - 2) 



for 
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Therefore, the regular solution of the Emden equation (|13|) behaves like the singular solution 
for £ — > +00. To determine the next order correction, we set z = zq + z' with z' <C 1. Keeping 
only terms that are linear in z', Eq. fllTD becomes 

(19) g + (D _ 2) £ +2( j>_ 2 y = . 

The discriminant associated with this equation is A = (D — 2)(D — 10). It exhibits two critical 
dimensions D = 2 and D = 10. For 2 < D < 10, we have 

^_2(D-2)/ 1 , A _J V /( J D-2)(10- J D) 



(20) e-* = v ^ y |l + p^cos^-* £ Mn£ + 5j j, (£ - +00). 

The density profile intersects the singular solution (JTJ) infinitely often at points that 
asymptotically increase geometrically in the ratio 1 : e 2 " :// V / (- D - 2 )( 1 °- D ) ( see; e .g. ; Fig. 1 of 
Ref. [0 for D = 3). For D > 10, we have 

.W, 2(D — 2) [ i 1 ^ ^ V(O-2)(10-£>) i -V(D-2)(10-D) 



(21) e-*= v ^ + s s ) >>, (< - vl 



For D = 2, Eq. (|T7| ) can be solved explicitly and we get 
(22) e"* 



" (1 + K 2 ) 2 ' 

This result has been found by various authors in different fields (see, e.g., {|, |22|). Note that 
~ £~ 4 at large distances instead of the usual £ -2 behavior obtained for D > 2. This implies 
that the mass of an unbounded isothermal sphere is finite in D = 2, although it is infinite for 
D > 2. 

For D < 2, we can neglect the r.h.s. in Eq. ( |T8D and we get 

(23) e-^~ e - A ^ 2 ~ D , (£^+00), 

where Ac is a constant depending on the dimension D. For D — 1, Eq. (|T3|) can be solved 
exactly, yielding the result 

(24) = ^ ^, 

V ' cosh 2 ^/^) 

establishing A\ = \/2. 

2.3 The Milne variables 

As is well-known [f2T |, isothermal spheres satisfy a homology theorem: if ^(0 is a solution 



of the Emden equation, then ip(A£) — 2 In A is also a solution, with A an arbitrary constant. 
This means that the profile of isothermal configurations is always the same (characterized 
intrinsically by the function ip), provided that the central density and the typical radius are 
rescaled appropriately. Because of this homology theorem, the second order differential equation 
(li~3| ) can be reduced to a first order differential equation for the Milne variables 



(25) 11 = — — and v = 
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> 2 




Figure 1: The solutions of the Emden equation in the (u, v) plane for systems with dimension 
2 < D < 10. 



Taking the logarithmic derivative of u and v with respect to £ and using Eq. (|13l) , we get 

1 du 1 



(26) 



(D-v-u), 



(27) 



1 dv 1 . _ . 

__ = _ (2 _ B+ll) . 



Taking the ratio of the foregoing equations, we obtain 

u dv 2 — D + u 
v du 



(28) 



D 



u 



The solution curve in the (u, v) plane is plotted in Fig. |I| for different values of D. The curve is 
parameterized by £. It starts from the point (u, v) = (D, 0) with a slope (dv/du) = — (D+2)/D 
corresponding to £ = 0. The points of horizontal tangent are determined by u — D — 2 and 
the points of vertical tangent by u + v = D. These two lines intersect at (u s , v s ) = (D — 2, 2), 
which corresponds to the singular solution fllip . For 2 < D < 10, the solution curve spirals 
indefinitely around the point (u s ,v s ). For D > 10, the curve reaches the point (u s ,v s ) without 
spiraling. For D = 2, we have the explicit solution v = 2(2 — u) so that (u,v) —>■ (0,4) for 
£ — > +oo. For D < 2, (u,v) — > (0, +oo) for £ — > +oo (see Fig. ||). More precisely, 



(29) 



ue< 



D 

V 2-D 



(£ - +oo) 



where we have defined uj d = 1/(A D (2 - D)^-o). For D = 1, ui = l/y/2. 



2.4 The thermodynamical parameters 



For bounded isothermal systems, the solution of Eq. ( |l3l) is terminated by the box at a normal- 
ized radius given by a = (SdPpoY^R- We shall now relate the parameter a to the temperature 
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Figure 2: The solutions of the Emden equation in the (u, v) plane for systems with dimension 
D = 1 and D = 2. 

and energy. According to the Poisson equation (01), we have 

(30) GM = J Gpd D v = j* S D Gpr D ~ l dr = j* j-(^ D - l ^jdr= (r^ 1 ^ ■ 
Introducing the dimensionless variables defined previously (using r/R = £/a), we get 

(31) V = ^D^I = a^{a). 

We note that, for D = 2, the parameter rj is independent on R. This is a consequence of the 
logarithmic form of the Newtonian potential in two dimensions. 

The computation of the energy is a little more intricate. First, using an extension of the 
potential tensor theory developed by Chandrasekhar (see, e.g., Ref. p3[), we can show that the 
potential energy in D-dimensions can be written 

(32) H/ = __i_ J pr .V$d D r, 

for D^2. Now, the Boltzmann-Poisson equation (|TTD is equivalent to the condition of hydro- 
static equilibrium 

(33) Wp = -pV$, 

with an equation of state p = pT. Substituting this relation in Eq. (|32|) and integrating by 
parts, we obtain 

(34) 2K+(D- 2)W = DV D R D p(R), 

where Vb = Sd/D is the volume of a hypersphere with unit radius. Eq. ([$4]) is the form of the 
Virial theorem in D-dimension. The total energy E = K + W can be written 

(35) E = + JL^V D R D p(R). 
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Figure 3: Evolution of the inverse temperature 77 along the series of equilibrium (parameterized 
by o) for 2 < D < 10. The curves correspond to D = 4, 3, 2.5, 2.2 from bottom to top. 

Expressing the pressure in terms of the Emden function, using p = pT and Eq. fll2"[) , and using 
Eq. ([H]) to eliminate the temperature, we finally obtain 

a = ER °~ 2 = D(*-D) 1 1 e^g 

1 j ~ GM 2 2(D-2)a^(«) £-2^'(o) 2 ' 

It turns out that the normalized temperature and the normalized energy can be expressed 
very simply in terms of the values of the Milne variables at the normalized box radius. Indeed, 
writing u = u(a) and v = v(a) and using Eqs. (pi]) (p6|), we get 

(37) 77 = v , 



(38) 



A 



1 

v 



~D(A-D) 
_2(D-2) 



Uq 



D 



The curves 77(a), A(a) are plotted in Figs. [||4]. For 2 < D < 10, they exhibit damped 
oscillations toward the values r) s — 2 and A s = 1/ (£) — 2) — £)/4, corresponding to the singular 
solution (0). For D > 10 the curves are monotonous. For D = 2, we have explicitly 



(39) 



?7 



or 



2(1 + |a 2 )' 



A 



a 1 



1 + 



o 



1 + 



ln 1 + 



o;- 



The expression of the energy has been obtained directly from Eq. (§) with the boundary condi- 
tion $(i?) = 0. The inverse temperature increases monotonically with o up to the value rj c = 4. 
Using Eq. fl22|) and returning to the original variables, we can write the density profile in the 
form 



(40) 



AM 



P 



7^(4-77X1 + ^ 
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In (a) 

Figure 4: Evolution of the energy A along the series of equilibria (parameterized by a) for 
2 < D < 10. 



This density profile is represented in Fig. |5| for different temperatures. At the critical inverse 
temperature r\ c = 4, all the particles are concentrated at the center of the domain. The density 
profile approaches the Dirac distribution 



(41) 



p(r) -> M5(r) for 



which has an infinite (negative) energy. 

For D < 2, the curves rj(a) and A(a) tend to +oo and as a 
explicitly 



+oo. For D — 1, we have 



(42) 



V 



2\/2atanh(a/v / 2), A 



1 



4v / 2atanh(a/ v / 2) 8 sinh 2 (a/V2) ' 



Note that according to Eq. (|32|), the energy is necessarily positive for D < 2, so the region 
A > is forbidden. Returning to original variables, the density profile is given by 



(43) 



M 



a. 



P 



tanh(a/V / 2) cosh 2 (^) 



where we recall that S% = 2. For a — > +oo, the profile tends to a Dirac peak MS(r). 

In Figs. 0-0, we have plotted the equilibrium phase diagram A — r], giving the temperature 
as a function of the energy, for different dimensions D. For 2 < D < 10, the curve spirals 
around the limit point (A s ,i] s ) corresponding to the singular solution. For D > 10, the curve 
is monotonous until the limit point. For D = 2, the curve is explicitly given by 



(44) 



A 



In 



and is represented in Fig. 0, together with the case D = 1. 

We stress that the preceding results, obtained in the meanfield approximation, are exact in 
the thermodynamical limit — > +oo such that rj and A are kept fixed. If the box radius is 
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Figure 5: Equilibrium density profile of a two-dimensional self-gravitating system as a function 
of the inverse temperature 77. For rj = 0, the density is uniform. For 77 — > r\ c — 4, the density 
tends to a Dirac peak. For r\ > r] c , there is no equilibrium state. 

given, this implies that T ~ N and E ~ N 2 . Alternatively, if the temperature and the energy 
per particle are given, the thermodynamical limit is such that N — > +00 with N/R D ~ 2 constant 
(for D > 2). 

2.5 The minimum temperature and minimum energy 

For 2 < D < 10, the curve 77(a) presents an extremum at points a n such that drj/da(a n ) = 0. 
Using Eqs. (|3T| ) and (0), we find that this condition is equivalent to 

(45) M = D - 2 = u s . 

Since the curve u = u s passes through the center of the spiral in the (u, v ) plane, this determines 
an infinity of solutions (see Fig. §), one a t each extremum of v (since r/ = t> ). Asymptotically, 



the a n follow a geometric progression (see Ref. |18| for more details): 

(46) a n ~ e 2 ™^^- 2 " 10 - ' (n -> +00, integer). 
In Fig. |3|, we see that an equilibrium state exists only for 

/ s PGM , s , 

(47) v = ^zj < iKai), (2 < ^ < 10). 

This determines a maximum mass (for given T and R) or a minimum temperature (for given 
M and i?) beyond which no equilibrium state is possible. In that case, the system is expected 
to undergo an isothermal collapse (see Sect. |3.3|) . For D = 2 and for D > 10, the 77(a) curve 
is monotonous. An equilibrium state exists only provided that 

(48) n = (3GM < Vc = 4, (D = 2), 
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Figure 6: Equilibrium phase diagram giving the inverse temperature rj as a function of minus 
the energy A for systems with dimension 2 < D < 10. 




Figure 7: Equilibrium phase diagram for two-dimensional self-gravitating systems. For infinitely 
negative energies, the inverse temperature tends to the value r\ c = 4. 
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Figure 8: Location of the turning points of energy and temperature in the (u, v) plane for 
systems with dimension 2 < D < 10. The construction is made explicitly for D = 3, which 
corresponds to the case extensively studied in Refs. fl7|, [Eg]. The dashed line v — 2 determines 



the location of the nodes of the density profiles that trigger the instabilities (see Sec. |2.6| ). 

/ s PGM 

(49) V = ^=2 <Vs = 2, (D > 10). 

We get comparable results for the energy. For 2 < D < 10, the curve A (at) presents an 
extremum at points a' n such that dA/da(a' n ) = 0. Using Eqs. fl3"8] ) and (|2T)|) Q2"7|) , we find that 
this condition is equivalent to 

(50) Au 2 + 2u v + (D 2 -8D + 4)u + D(D - 2) (4 - D) = 0. 

We can check that the limit point (w s , w s ) is solution of this equation. Therefore, the intersection 
of the parabola (P) defined by Eq. ([5(]) with the spiral in the (u, v) plane determines an infinity 
of points a' n at which the energy is extremum (see Fig. |8]). On Fig. f|, we see that an equilibrium 
state exists only for 

-ER D - 2 

(51) A= — <AK), (2< J D<10). 

This determines a minimum energy (for given M and R) or a maximum radius (for given M and 
E) beyond which no equilibrium state exists. In that case, the system is expected to collapse 
and overheat; this is called gravothermal catastrophe (see Sect. |3.4j ). For D > 10, the curve 
A(a-) is monotonous. An equilibrium state exist only for 

rp tdD—2 i tj 

(52) A = Z_< A , = ___ (C > 10 ). 

For D = 2, there exists an equilibrium state for each value of energy (see Fig. 0): there is no 
gravothermal catastrophe in the microcanonical ensemble in two- dimensions Q. For D < 2, 
there exists an equilibrium state for all accessible values of energy and temperature. 
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2.6 The thermodynamical stability 



We now study the thermodynamical stability of self-gravitating systems in various dimensions. 
We start by the canonical ensemble which is simpler in a first approach. A critical point of 
free energy at fixed mass and temperature is a local maximum if, and only if, the second order 
variations 



(53) 



6 2 J 



2p 2T 



5p5M D r 



are negative for any variation 5p that conserves mass to first order, i.e. 



(54) 



Spd D v = 0. 



This is the condition of thermodynamical stability in the canonical ensemble. Introducing the 
function q(r) by the relation 



(55) 



6p 



dq 



SdT d 1 dr ' 



and following a procedure similar to the one adopted in Ref. ||18|1 , we can put the second order 
variations of free energy in the quadratic form 



(56) 



6 2 J=~ 



R 



drq 



G 



+ 



d ( 



d 



Tr D 1 dr\Srjpr D 1 dr 



q- 



The second order variations of free energy can be positive (implying instability) only if the 
differential operator which occurs in the integral has positive eigenvalues. We need therefore 
to consider the eigenvalue problem 



(57) 



d 



d 



drySopr 1 dr J Tr D 



G 



A<?a(?") 



with q\(0) = q\(R) = 0. If all the eigenvalues A are negative, then the critical point is a 
maximum of free energy. If at least one eigenvalue is positive, the critical point is an unstable 
saddle point. The point of marginal stability, i.e. the value of a in the series of equilibria 
r](a) at which the solutions pass from local maxima of free energy to unstable saddle points, is 
determined by the condition that the largest eigenvalue is equal to zero (A = 0). We thus have 
to solve the differential equation 



(58) 



d 



dF 



dr \SDpr D 1 dr J Tr D 



+ 



GF 



0. 



with F(0) = F(R) = 0. Introducing the dimensionless variables defined previously, we can 
rewrite this equation in the form 



(59) 

with F(0) = F(a) 
(60) 



1D-1 



0. 



0. If 



d ( d \ 



-.D-l 
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denotes the differential operator which occurs in Eq. (|59|), we can check by using the Emden 
Eq. (0) that 

(61) C(e~ V) = £(£ D e^) = (D - 2W. 

Therefore, the general solution of Eq. (|59D satisfying the boundary conditions at £ = is 

(62) F(0 = Cl (£ D e^ - (D - 2)^V). 

Using Eq. (|62|) and introducing the Milne variables (^51), the condition F(a) = can be written 

(63) uq = D — 2. 

This relation determines the points at which a new eigenvalue becomes positive (crossing the 
line A = 0). Comparing with Eq. (f45l), we see that a mode of stability is lost each time that rj is 
extremum in the series of equilibria, in agreement with the turning point criterion of Katz |20| 
in the canonical ensemble. In particular, the series becomes unstable at the point of minimum 
temperature (or maximum mass) a±. Secondary modes of instability appear at values «2, 
We obtain the same results by considering the dynamical stability of isothermal gaseous spheres 
with respect to the Navier-Stokes equations (see Ref. [TS|] for D — 3). Therefore, dynamical 
and thermodynamical stability criteria coincide for isothermal gaseous spheres. 

According to Eq. ([55]), the perturbation profile that triggers a mode of instability at the 
critical point A = is given by 

(64) S ± = -J—'^l 
{ ' Po Sue- 1 di ' 

where F(£) is given by Eq. fl6"2"|). Introducing the Milne variables (|25|), we get 

(65) ^ = f(2-«). 

The density perturbation Sp becomes zero at point(s) such that f(£j) = 2. The number of 
zeros is therefore given by the number of intersections between the spiral in the (u, v) plane and 
the line v = 2 (see Fig. [8]). For the n-th mode of instability we need to follow the spiral to the 
n-th extremum of v (since a n corresponds to an extremum of rj, hence v). Therefore, the density 
perturbation 5p corresponding to the n-th mode of instability has n zeros Ch^2, ■■■,Cn < ««■ 
Asymptotically, the zeros follow a geometric progression with ratio 

In the microcanonical ensemble, the condition of thermodynamical stability requires that 
the equilibrium state is an entropy maximum at fixed mass and energy. This condition can be 
written 

(66) ?s = _y <|V r _ _l / Wr _ _L_(y w ^ r y < o, 

for any variation bp that conserves mass to first order (the conservation of energy has already 
been taken into account in obtaining Eq. (|66|)). Now, following a procedure similar to that of 
Ref. [§] in D = 3, the second variations of entropy can be put in a quadratic form 

(67) 8 2 S = / / drdr'q(r)K(r,r')q(r'), 



o 
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with 

< 68 > K ^ = -D^ d ^§^ +1 2^-^ 



G d ( 1 



Tr D ~ l dr\SDpr D ~ 1 dr / 
The problem of stability can therefore be reduced to the study of the eigenvalue equation 



(69) 



I dr'K{r,r')F x {r') = XF x {r) 
Jo 



with -Fa(O) = F X (R) = 0. The point of marginal stability (A = 0) will be determined by solving 
the differential equation 

d ( 1 dF\ GF 2V d^ 
( } d?\S D pr D - l lr~ J + T^ 1 ~ DMT 2 ~dr 

with 

(71) V= / — {r')F{r')dr' . 

Jo dr 

Introducing the dimensionless variables defined previously, this is equivalent to 

d ( dF\ F _ dtp 



with 

and -F(O) = F(a) = 0. Using the identities (|6T|), we can check that the general solution of 
Eq. (fr^) satisfying the boundary conditions for £ = and £ = a is 

(74) F(0 = 7^ (Vv* - (£> - 2)£ D ~ y) + X^" V, 

The point of marginal stability is then obtained by substituting the solution ( f74[) in Eq. (|73|). 
Using the identities 

(75) f rf/^e^ = V (a)(£> - u ), 



(76) (£>-2) /V~W^ = ^~V(a)(2£>-2tzo-Vo), 

./o 

which result from simple integrations by parts and from the properties of the Emden equation 
(fl3|), it is found that the point of marginal stability is determined by the condition (|50|). There- 
fore, the series of equilibria becomes unstable at the point of minimum energy in agreement 
with the turning point criterion of Katz [f2D| in the microcanonical ensemble. The structure of 
the perturbation that triggers the instability can be determined with the graphical construction 
described in Ref. ||17| . It is found that the first mode of instability has a core-halo structure 
(i.e., two nodes) in the microcanonical ensemble, unlike the first mode of instability in the 
canonical ensemble JTS[ . 

The thermodynamical stability analysis presented in this section also shows that the equi- 
librium states for D < 2 and D > 10 are always stable since the series of equilibria do not 
present turning points of energy or temperature. 
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3 Dynamics of self-gravitating Brownian particles in di- 
mension D 



3.1 The Smoluchowski-Poisson system 

We now consider the dynamics of a system of self-gravitating Brownian particles in a space of 
dimension D. Like in Ref. Q, we consider a high friction limit in order to simplify the problem. 
We thus study the Smoluchowski equation 



(77) g = V 



i(TVp + pV$) 



coupled to the Newton- Poisson equation (Q). In the microcanonical ensemble, the temperature 
T(t) evolves in time so as to satisfy the energy constraint (|]). In the canonical ensemble, 
the temperature is constant. The Smoluchowski equation can be obtained from a variational 
principle called the Maximum Entropy Production Principle 0. This variational approach is 
interesting as it makes a direct link between the dynamics and the thermodynamics. In the 
microcanonical description, the rate of entropy production can be put in the form |2[] 

(78) S= [ -^—(TVp + pV$) 2 d D r > 0. 



For a stationary solution, S = and we obtain the Boltzmann distribution fllOD which is a 
critical point of entropy. Considering a small perturbation around equilibrium, we can establish 
the identity §: 

(79) 5 2 S = 2X5 2 S > 0, 

where A is the growth rate of the perturbation defined such that 5p ~ e xt . This relation shows 
that a stationary solution of the Smoluchowski-Poisson system is dynamically stable if and only 
if it is a local entropy maximum. We get similar results in the canonical ensemble with J in 



place of S. The relation ( |79|) has been found for other kinetic equations satisfying a if-theorem 
(Chavanis, in preparation). In D = 2, an equation morphologically similar to the Smoluchowski 
equation (|77|) can been derived from the iV-body Liouville equation for a gas of point vortices 
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3.2 Self-similar solutions of the Smoluchowski-Poisson system 

From now on, we set M — R — G — £ — 1. The equations of the problem become 

(80) ^ = V(TVp + P V$), 

(81) A$ = S D p, 

(82) e=^T+ 1 -! p$rf D r, 



2 2 



with boundary conditions 



9$ 1 dp 

(83) ^(0,t) = 0, ®( 1 ) = ^Td> T ^( 1 ) + ^( 1 ) = ' 
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for D > 2. For D = 2, we take $(1) = on the boundary. Integrating Eq. fl3T|) once, we can 
rewrite the Smoluchowski-Poisson system in the form of a single integrodifferential equation 

The Smoluchowski-Poisson system is also equivalent to a single differential equation 
, x dM rT1 fd 2 M D-ldM\ 1 fiM 

(85) -5T = T U— - — — + —jjr—rrM- 



for the quantity 



(86) M(r,t) 



/ p(r')SDr D ~ 1 dr', 
Jo 



which represents the mass contained within the sphere of radius r. The appropriate boundary 
conditions are 

(87) M(0,t)=0, Af(l,t) = l. 

It is also convenient to introduce the function s(r,t) = M(r,t)/r D satisfying 



ds _ T fd 2 s D + l ds\ fds \ 
dt \dr 2 r dr J \ dr J 



We look for self-similar solutions of the form 



T \ 1/2 



(89) p(r, t) = po(t)f I j— | , r = - 

VoWJ \Po 

In terms of the mass profile, we have 

(90) M(r,t) = M (t)g(^-^J, with M {t) = p r D , 
and 

(91) g(x) = [ f(x')S D x' D - l dx'. 

Jo 

In terms of the function s, we have 



(92) s(r , t) = po(t)S I — t-t I , with S(xj- n 

\?WJ/ x 

Substituting the ansatz d92|) into Eq. ([S8|), we find that 

(93) ^S(x) - ^^x5'(x) = pi (s"{x) + ^V(x) + xS(x)S'(x) + DS(x 



where we have set x = r/r . The variables of position and time separate provided that there 
exists a such that 

(94) p r« = K, 
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where k is a constant. In that case, Eq. (|93|) reduces to 

(95) ^ (s(x) + -xS'(x)) = pi (s"(x) + E+ls'(x) + xS(x)S'{x) + DS(x) 2 

Assuming that such a scaling exists implies that (1/ pfy^dpo/dt) is a constant that we arbitrarily 
set equal to a (note that this convention is different from the one adopted in Ref. This 
leads to 

(96) Po (t) = -{t coll -t)-\ 

a 

so that the central density becomes infinite in a finite time t co u. The scaling equation now reads 

(97) aS + xS' = S" + ^-^S' + S{xS' + DS). 
For x — > +oo, we have asymptotically 

(98) S(x) ~ </(x) ~ x D " a , /(a;) ~ x~ a . 
3.3 Canonical ensemble 

In the canonical ensemble in which the temperature T is a constant, we have [] 

(99) a = 2, k = T. 



In that case, the scaling equation (|97|) can be solved analytically Following a procedure similar 
to the one developed in Ref. [0], we find that 

A 

(100) S{x) 



D-2 + x 2 
Then, Eqs. (|92]) and (fH|) yield 



(101) (/(a;) = - - ■ „ /(z 



4(D-2) L> + x 2 



D-2 + x 2 ' JK >~ S D ( jD _ 2 + x2) 2 ' 

According to Eqs. (|89| ) and (0), the central density evolves with time like 

2D 

(102) P (o, t) = p (t)/(o) = _ (t coll - ty 1 . 

According to Eqs. ( [89]) and (0), the core radius and the core mass evolve like 

(103) r (f) = v^t**, - t) 1 / 2 , M„(f) = i(2T)^ 2 (t co « - f)*" 1 . 

Note that for D > 2, the core mass goes to zero at the collapse time. At t — t co u, we get the 
singular profile 

(104) p(r, t = t coll ) = AT ^~ 2 \ M(r, t = t coll ) = 4Tr D ~ 2 . 

S D r 2 



1 The case T = is treated in the Appendix 
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3.4 Microcanonical ensemble 



In the microcanonical ensemble, the exponent a is not determined by simple dimensional anal- 
ysis. In Ref. JIJ, we found numerically that the scaling equation ( p7|) has physical solutions 
only for a < a max , with a max ~ 2.21 for D = 3. We also argued that the system will select the 
exponent a max , since it leads to a maximum increase of entropy. In this section, we show that 
in the limit of large dimension D, we can explicitly understand the occurrence of such a ct max . 
In addition, we will present the derivation of perturbative expansions for ct max and the scaling 
function S, in powers of D~ L . 

Eq. ( p7|) can be formally integrated as a first order differential equation (writing S" = 
S' x [S"/S']), leading to an expression of S(x) as a function of x, S(x) itself, and S"(x)/ S'(x): 



(105) 



a 



DS(x) 



- 1 



a 



DS(0) 



- 1 



exp 



a 



ydy 



o y 2 (l-S(y))-y§§ 



(D + l) 



We now define x , such that S(x ) = a/D. Since S should be analytic, the foregoing relation 
implies for x — > xq, 



(106) 



ay 
o Hv) 



dy ~ In \x — Xq\ 



where F(y) is the function which occurs in the denominator of the integral in Eq. ( |105| ). From 
Eq. ( |106|) , we must have F(y) = axo(y — Xq) for y — > Xq, which implies F(xq) = and 
F'(xq) = axQ. These conditions can be rewritten 



(107) 



0, 



108) 



[a 



2)x = -± 
ax 



x 2 S(x) + x 



S'(x) _ 



This preparatory work now allows the introduction of a systematic expansion in large dimension 
D for the scaling function S, the scaling exponent a, and xq. In this limit, let us neglect the 
contribution of the terms which are not of order D in the right-hand side of Eq. (0). This 
actually amounts to take F(y) ~ y 2 — D in Eq. ( |105| ). Within this approximation, we find 



(109) 



a 



DS(x) 



a 



1 



x 



x 
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which is an analytic function only if a 
form for S, 



(110) 



S(x) 



DS(0) 
■ 2. This leads to x 

3(0) 
l+(«-l 



a/2 



D, and to the more explicit 



v 



S(0) remains undetermined, and will be fixed by the next order approximation. Indeed, we can 
iteratively solve the full scaling equation Eq. (pTj), by reinserting the zeroth order solution into 
Eq. (|105|) , and eventually continue this process with the new improved scaling function, and 
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so forth... Thus, expressing the conditions of Eq. (|107| ) and Eq. (|108|) , and defining z 
(which will be of order 0(1)), we obtain 



DS(0) 

2 



(111) 

and 
(112) 



x 2 = D + - + O^D- 1 ), O tx = Vd(i + ^- + 0(D- 2 )] 
z \ zD J 



a 



4 

2 = — 
D 



z z^ 



0{D~ 



Eq. ( |112[ ) provides a relation between the possible values for a and the associated value of 
S(0) = j}. Note that the function of z in the right-hand side of Eq. (|112|) has a well defined 
maximum. Hence, up to order 0(D^ 1 ), we find that a < a max , with 



(113) 



a m&x = 2+ l -D- 1 + 0{D- 2 ) 



which is associated to the value z = 4 + 0(D _1 ) or S(0) = 4 + 0(D~ 2 ). As a is necessarily 
greater than 2 (as the temperature cannot vanish), a solution exists for any a G [2, a max ]. As 
already mentioned, a max is dynamically selected as it leads to the maximum divergence of the 
entropy and the temperature (see Eq. (|123|) below). 

Inserting Eq. ( |110| ) into Eq. (|105|) , we find the next order approximation for S 



(114) 



a 



DS(x) 



a 



DS(0) 



x 



x 



.if 



f(i-« 



X, 



+ 1 



where xq is given by Eq. (|111|) , and x\ and are defined by 



(115) 



D 



z{z-l) 



2 
D 



z z' 



+ 0{D~ 



Again, the analyticity condition imposes that ^(1 
and to the following explicit form for S: 



1, which exactly leads to Eq. (|112 ), 



(116) 



S{x) 



a 
D 



1+1 



a 



X 



.If. 



X 



.) I 



r-1 



This improved scaling function can be inserted again into the conditions expressed by Eq. ( |107| ) 
and Eq. (|108|) , leading to the next order term in the expansion of a. After elementary but 
cumbersome calculations, we end up with 



(117) 

a - 



2 = 1 
D 



+ D 2 



5 26 31 6 

z z 2 z 3 



z z' 

This function has again a well defined maximum for 



1 7 14 8 . , 
^ + ^r r I m z 

z Z A Z A z 



+ 0{D 



-3\ 



(118) z 

associated to the value 

(119) a t 



jS(0) = 4 + - 61n2^) D- 1 + 0(ZT 2 ), 



2 + -D- 1 + —D- 2 + 0(D- 3 ). 
2 16 K ' 
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This expansion gives a max — 2.24... in D = 3, in fair agreement with the exact value a max = 
2.2097... obtained numerically in ][]]. In addition, the exponent a = 2 is associated to z — 
2 + 4_D~ 1 + 0(D~ 2 ). In principle, these expansions can be systematically pursued to the prize 
of increasingly complicated calculations. 

Finally note that Eqs. fl39] ) and ([)6|) lead to the following exact asymptotic for the central 
density p(0, t): 

2z(a) 

(120) p(0,t)~K D (a)(t coll -t)-\ K D {a) 



aS 



D 



where we have used /(0) = DS(0)/Sd and the definition of z. The function z[a) is determined 
implicitly by Eq. ( |117|) , up to order 0(D~ 2 ). For the special cases a = 2 and a = a max , we 
respectively find 

(121) K D {2) = 2S^ 1 (1 + 2ZT 1 + 0(ZT 2 )) , 

(122) #j>(<w) = 45 D 1 (l+^-^ln2^- 1 + 0( J D- 2 )^), 

which shows that i^D(a ma x) is substantially greater than Kd(2) (twice bigger in the infinite D 
limit, the ratio being even bigger for finite D, as y — fin 2 « 3.835... > 2). This substantial 
difference was noted in Ref. |J, in the case D = 3. Finally, as expected in the microcanonical 
ensemble, the temperature diverges during the collapse as T(t) ~ {t co ii—t)~ aT with «t = 1— 2/a; 



see Eqs. (|39"|) (|94|) (p6|) . The strongest divergence is obtained for a = a max . According to 



Eq. ( |119| ), we have 



(123) a T = 2 = t -D + — D~ + 0(D 

a max 4 32 

If we plug D = 3 in Eq. ( p. 231 ), we fi n d the estimate «t w 0.11... in fair agreement with the 
exponent measured numerically in olt ~ 0.1. 

4 The two-dimensional case 
4.1 The critical temperature 

In two dimensions, the dynamical equation ( P5| ) for the mass profile reads 

(124) ir = 4r "^ + 2M *- 

after the change of variable u = r 2 has been effected. Looking for a stationary solution, and 
using uM" = (uM'Y — M', Eq. (|124|) is readily integrated leading to 



AT u 

(125) M{u) 



AT -I 1 + 



Note that M(l) = 1, ensuring that the whole mass is included in this solution. Using p = M' /n, 
we find that the density profile is given by 

(126) P (r) AP ° 1 



vr (l + ( r /r ) 2 ) 5 
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with 



(127) r = V4T-1 and p r 2 = T. 

This solution exists provided that T > T c = 1/4, which defines the collapse temperature. We 
have thus recovered the result ( fHf ) by a slightly different method. Note that the value of T c 
and the dependance of Tq and po with the temperature coincide with the exact results obtained 
within conformal field theory In the following, T is set constant as we have already seen 
that the gravothermal does not exist in the microcanonical ensemble in two dimensions. 



4.2 Scaling collapse for T = T c 

We now address the dynamics at the critical temperature T = T c = 1/4. We note that contrary 
to what happens in other dimensions, the central density diverges at T c . Thus, in analogy with 
critical phenomena, we anticipate a scaling solution for M(u,t), of the form 

(ct(t) + l)u 

128 M(u, t) J , 

1 + a(t)u 

which preserves the scaling form obtained above T c , and which satisfies the boundary condition 
M(l, t) = 1. The corresponding density profile is 

/ x a(t) + 1 1 
(129) p(r,t) = 



7T (l + a(t)r 2 Y 
The central density 

a(t) + 1 



(130) p(0,t) 



71 



is expected to diverge for t — > +oo, so that a(t) is also expected to diverge. 

Inserting the ansatz Eq. ( |128|) into Eq. ( |124| ) shows that the left-hand term is indeed negligi- 
ble compared to both terms of the right-hand side, to leading order in a. So far, this prevents us 
from determining a dynamical equation for a. In order to achieve that, we must solve Eq. ( |1 24j ) 
to the next order in a -1 . We thus look for a solution of the form 

(131) M(u, t) = + a{ty l h{u, t), 

1 + a(t)u 

where h(u, t) is expected to be of order 0(1), and should satisfy h(0, t) = and h(l, t) — 1 (the 
total integrated mass should be and 1, respectively, for u = and u = 1), and |^(0, t) = 0, 



which ensures that Eq. ( |130| ) is exactly obeyed, defining a(t) without any ambiguity. The 
contribution of in the left-hand side of Eq. ( 124 ) is dominated by the time derivative of 
Eq. O : 



< 132 > jft 



[1 + a(t))u 
1 + a(t)u 



u(l — u) da 
(l + au) 2 ~dl'' 



which will be checked self-consistently hereafter. In addition, non linear terms in h in the 
right-hand side are also negligible. Therefore, h satisfies 



au 



au(l — u) da d 2 h d , 
133 — i -f— = u— + 2— h 

K J 1 + au 2 dt du 2 du \ 1 + au 
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Figure 9: At T — T c = 1/4, and when the central density has reached the value p(0,t) ~ 
1644.8... = a ^ +1 (a(i) ~ 5166.3...), we have plotted the result of the numerical calculation 
compared to our exact scaling form p(r,t) = a ^ +1 (1 + a{t)r 2 )~ 2 obtained in Eq. (|129| ). The 



two curves are indistinguishable as the relative error is, as predicted, of order a ~ 10 . Note 
finally that for this range of density, the density contrast is huge, of order 10 7 . 



Actually, for a given time, this equation becomes an ordinary differential equation involving 
only one variable u, as a and da/dt appear as parameters. Eq. (133) can be integrated leading 
to a first order equation in h, which can be easily solved. Defining v = au, we finally get 



(134) 



h(u, t) 



2\ da 



a ) dt 



:i + v) 



X 



l)ln(l + v) +v(l - 2v) +2v 

Jo 



v ln(l + z) 



dz 



2v 2 + v 3 
2(a + 2) 



which depends on time only through the variable a and da/dt. Now, da/dt is determined by 
imposing the boundary condition h(l,t) = 1, which leads to 



(135) 



da 
~dl 



In a - 5/2 



[l + 0(lna" 2 )] . 



One can solve iteratively Eq. ( 124|) , by adding the time derivative of the above solution to the 
left-hand side, in order to compute an improved h. To leading order, the solution of Eq. ( |134| ) 
is preserved. However, new terms are generated which are important for v ~ a (u ~ 1), and 
which generate terms of order 0(a/ In 3 a) is the expansion for da/dt. This explains the form of 
the error term in Eq. (|135|) . 

Integrating Eq. ( |135| ) for large time, we get the exact asymptotic expansion for large time 



(136) 



a{t) = exp Q + V2tj [1 +0(r 1/2 lnt)] . 



For t — ► +oo, the central density diverges like a(t) and the core radius goes to zero like a(t) 1 ^ 2 . 
In addition, the scaling solution ( |129|) at T = T c goes to a Dirac peak containing the whole 
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2 4 6 8 

ln(a) 

Figure 10: We plot alda/dt) -1 as a function of In a, which is predicted to behave as 
a(da/dt)~ l ~ In a — 5/2 + 0([lna] _1 ) (see Eq. ( |135| )). Even for the moderate range of ac- 
cessible densities (a max ~ 5166), we clearly find that the numerical result evolves toward the 
theoretical asymptotic (dashed line). 

mass (see Eq. (pH])), as the decay exponent of the scaling function is 4, which is strictly greater 
than 2. 



4.3 Collapse for T <T C 



For D = 2, the scaling equation associated to Eq. (|85|) does not display any physical solution 
when numerically solved. In this section, we thus present a special treatment adapted to this 
case. The principal difference with other dimensions is the divergence of the central density at 
T c , and the occurrence of a scaling solution at this temperature. 

Strictly below T c , we expect a finite time collapse. Close to the center, the solution takes 
the form 

a(t)u 

(137) M(M)~4T : 



l + a(t)u 

where again the left-hand side of Eq. ( |124j ) is negligible compared to each term on the right-hand 
side. We thus look for a solution of the type 

(138) M(u, t) = AT a } t] l + h(u, t), 

1 + a(t)u 

where h is now of order 0(1) as it contains a finite fraction of the total mass, since the first term 
contains a mass of order 4T < 1. Inserting this ansatz in the dynamical equation Eq. ( |124J ), we 
obtain 

1 dh da u d 2 h d ( au \ dh 

(139) 4T dt + ~dt(l + auY = U drf + 2 frl [ T+a^t + 
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One can look for a scaling solution of the type 

(140) h(u,t) = a i ~ x H(au), with H(v) ~ cv 1 ^, when v 



-oo, 



so that the mass included in this scaling profile h(l,t) = c = 0(1). With this definition, the 



density profile decays for large distance as p 
Eq. ( |139| ), we obtain 



with a = 27. Inserting this ansatz in 



(141) 



^ (^ + (7 - 1)# ) W^-^ 



—a = vH 
dt 



where derivatives are with respect to the variable v. We are free to choose a(t) = np(0, t)/(4T), 
so that H'(0) = H(0) = 0. For small v, Eq. ([TIlD leads to 



(142) 



^ = H"(0)a^ +1 . 
dt y ' 



Eq. ( |141| ) has a global scaling solution only for 7 = 1. However, we know that in this case the 
scaling equation obtained by setting 7 = 1 does not display any physical solution. Thus, we 
conclude that there is no scaling solution obtained by imposing that all terms in Eq. (|141|) scale 
the same way. However, as we will see in the section devoted to numerical simulations, the direct 
simulation of Eq. ( |124| ) seems to display a scaling solution with 7 « 0.6 — 0.7 for numerically 
accessible densities. Strictly speaking, this is totally excluded by the above equation, except 
if one allows 7 to very slowly depend on the density or a. For a given a, we thus want to 
solve Eq. (|141|) , where the boundary conditions will ultimately select the effective value of 7, 
and that of da/dt. More precisely, once we impose H'(0) = H(0) = 0, and the condition of 
Eq. (|142| ), we end up with a shooting problem for H"(0) and 7. For large a, and v <C a, it is 
clear that the non linear term of the right-hand side of Eq. ( |141| ) becomes irrelevant, and we 
drop it from now on. 

In order to understand the origin of this shooting problem, and to obtain an accurate 
estimate of 7, let us solve Eq. ( |141| ) in the limit of very large a, in the range 1 <C ?; <C a. In 
this regime, Eq. (|141|) simplifies to the following equation 



(143) 

where 
(144) 



1 

4T 



vH' + (7 - l)H) + a 1 -^- 1 



CO 



vH" + 2H', 



UJ 



da 
~dt ( 



H"(0)a 



7-1 



Let us now multiply this equation by v 7 2 and integrate the resulting equation. After elemen- 
tary manipulations, we obtain 



(145) 
H' + 



3 — 7 00 
~v 4T 



H = -^v^ - ^V 1 + (2 - 7 )(3 - 7)^ 
4i 2 — 7 



w J ~ 3 H(w) dw. 



where c ~ 0(1), which has been defined in Eq. ( |140|) , appears here as an integration constant. 
Then, one can integrate this differential equation which leads to the following self-consistent 
relation for H: 



(146) 



/-|-oo 
■u; 3 ~ 7 exp 



oow\ 

■—— F(w) dw, 
AT J v 1 ' 
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where F is defined as the opposite of the right-hand side of Eq. ( |145| ) : 



toe 1 „, uja 1 ^ r+™ 



(147) F(v) = ^v 1 -i+ v~ l - (2- 7 )(3-7)^ 1 ~ 7 / w^ 3 H(w)dw. 

4T 2 - 7 A 

Eq. ( |143|) implies that -ff(f) ~ lnt>, when v — > (of course, this apparent divergence does not 
occur in the full dynamical equation Eq. (|141|) ). Considering the prefactor t> 7-3 in Eq. (|146|) , 
this behavior can be obtained if and only if 

(148) l + w 3 - 7 exp(-^) F(w)dw = 0. 

J 



As u is expected to go to zero for large a as 7 < 1 (see Eq. (144)), the dominant contribution 
of the integral of the third term in the definition of F comes from the large w region, for which 
H can be replaced by its asymptotic form. Hence, defining T(x) = / +o ° w x exp(-w) dw and 



e = 1 — 7, and using Eq. ([[44]), the condition expressed in Eq. ( |148|) can be rewritten 



(149) C = „ 



e(l + e) 2 T(l + 2e) w V 4T 



e 



As c is of order 0(1), we find that e — > as a — > +00. More precisely, in this limit, e is the 
solution of the following implicit equation 



(150) 



\n(K/e) 



In a 

where K = H"(0)/c + 0(e). Finally, we obtain 



(151) e = 1 - 7 = (1 + Oainlna]- 1 )) . 

In conclusion, although the solution is not strictly speaking a true scaling solution, the explicit 
dependence of 7 on a is so weak that an apparent scaling should be seen with an effective 7 
almost constant for a wide range of values of a. Hence, the total density profile is the sum of 
the scaling profile obtained at T c with a T/T c weight (behaving as a Dirac peak of weight T/T c , 
at t = t co u) and of a pseudo-scaling solution associated to an effective scaling exponent slowly 
converging to a = 2. 

Let us illustrate quantitatively the time dependence of a = 27. For example, taking arbitrar- 
ily K = 1 (the dependence on K is weak and vanishes for large a) , Eq. ( |150| ) and Eq. ( |151| ) re- 
spectively lead to 7(0 = 10 3 ) = 0.624... and 7(0 = 10 3 ) = 0.626..., and to 7(0 = 10 5 ) = 0.684... 
and 7(0 = 10 5 ) = 0.674... (note that the error between the asymptotic expansion of Eq. ( |151| ) 



and the implicit expression first grows before slowly decaying for a 3> 10 12 !). Finally, for the 
maximum value of a accessible numerically of order a ~ 10 4 , we expect to observe an apparent 
scaling solution with 7 0.65, or a = 27 w 1.3. 



4.4 Numerical simulations 

In this subsection, we present numerical simulations concerning the two-dimensional case. In- 
deed, the three-dimensional case has been extensively studied in [0. It has been shown that 
the scaling function as well as the corrections to scaling (which have been calculated for the 
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10" 2 10" 1 10° 10 1 10 2 

X=r ^ r max 



Figure 11: At T = T c /2 = 1/8, we have extracted the next correction to scaling p cor = p — 
4TpT=T c , where Pt=t c is defined in Eq. ( |129| ). We have then plotted p cor {r, t) / p cor (r max (t) ,t) as 
a function of x = r/r max (t), where r max (t) is defined as the location of the maximum of p CO r(r, t). 
Consistently with the apparent scaling observed, we found r^l x (t) ~ *Ja ~ a/ p C or( r max(^), t). 
For a = 2 n ~ 1 x 100 (n = 1,...,8), we have obtained a convincing data collapse associated to 
a — 2j « 1.3, in agreement with the theoretical estimate of 7, in this range of a. 



canonical ensemble in U) are perfectly described by the theory In addition, as the system 
behaves qualitatively the same for any dimension D > 2 in both thermodynamical ensembles, 
we have decided to focus on the numerical study of the D = 2 case, which displays some very 
rich behaviors, as exemplified in the present section. 

We consider the system in the canonical ensemble, as the gravitational collapse does not 
occur in the microcanonical ensemble. In Fig. |9|, we show the scaling function at T c , as given 
by Eq. ( |129| ), finding a perfect agreement. In Fig. [TU], we also display a(da/dt)~ l as a function 
of In a, and find an asymptotic behavior fully compatible with that given by Eq. ( 135|) . 

Below T c , and in the accessible range of a (up to a ~ 10 5 ), we find an apparent scaling regime 



with a = 27 m 1.3, as predicted in Sec. [0| This is illustrated in Fig. |TT], for T = T c /2 = 1/8. 
Note that the effective 7 or a can also be extracted from the time evolution of a(t) or the central 
density (see Eq. ( |142| )). In Fig. [12], we show that this way of measuring 7 is fully compatible 



with the value of the effective exponents a — 2j « 1.3 



5 The one dimensional case 

When an equilibrium state exists, there is little hope to be able to solve the full Smoluchowski- 
Poisson system analytically in order to study the relaxation towards equilibrium. We shall 
consider a simpler problem in which a test particle evolves in a medium of field particles at 
statistical equilibrium (thermal bath approximation). The particles are assumed to create a 
stationary potential $eq(r) which induces a drift of the test particle along the gradient of $ e(? . 
In addition, the test particle is assumed to experience a diffusion process. If p denotes the 
density probability of finding the test particle in r at time t, we expect the evolution of p to be 
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Figure 12: We plot aT l (da/dt) ~ a 1 as a function of a, in order to extract the effective value of 
7 directly from the time evolution of the central density. We find that the effective 7 is slowly 
growing with time, as predicted, and is of order 7 = a/2 pa 0.65 (the dashed line has a slope 
equal to 0.65), which is fully compatible with the value extracted from Fig. D], and the value 
expected from Eq. (|150|) in this range of a. 

determined by a Smoluchowski equation of the form 



(152) ^ = V(TVp + pV$ e ,), 



where $ e <?( r ) is solution of the Boltzmann-Poisson equation (|ll|). This means that we replace 
the true potential by its equilibrium value but still allow the density p to vary with time. As 
we shall see, it is possible to solve the Smoluchowski equation ( |152j ) analytically in D = 1 by 
using an analogy with a problem of quantum mechanics. An equation of the form (|152 ) has 



been proposed in Refs. [p7| , |25|j to model the motion of a test vortex in a bath of field vortices 
at statistical equilibrium. In that context, Eq. (|152|) can be formally derived from the A-body 
Liouville equation of the system by using projection operator technics. 

It is well-known that a Fokker-Planck equation like Eq. ( |152j ) can be formally transformed 
into a Schrodinger equation with imaginary time. Indeed, performing the change of variable 

(153) p = ^ e -^* e % 

we find that the evolution of ip is determined by an equation of the form 

(154) ^=TA^+ Qa$ 69 - ^(V$ e<? ) 2 ^. 
This can be written as a Schrodinger-type equation 

(155) ^ = TAip — V(r)i>, 
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with a potential V(r) = — ^A& eq + ^(V$ e g) 2 . So far, this transformation is general. If we now 
consider the one-dimensional case, the Boltzmann-Poisson equation ( [L52| ) can be solved ana- 
lytically and the potential V(r) determined explicitly. Introducing the notations £ = ar/^/2R 
and r = a 2 Tt/2R 2 and using Eq. (El) we can rewrite Eq. (154) in the form 



< 156 > f^ + fe- 1 )*' 

A separation of the variables can be effected by the substitution 

(157) ^(e,t)=0(Oe" A ' (A>0), 
where <fi is solution of the ordinary differential equation 

(158) V + 2[E + 



d£ 2 \ cosh' £ 

where we have set A — 1 = 2E. The solutions of this Schrodinger equation are described in 
detail in Ref. ||28|| . The spectrum of positive energies is continuous. The spectrum of negative 
energies is discrete and reduces to Eq = —1/2 (fundamental state). The first excited state in the 
continuum is E\ = 0. We can check that the corresponding eigenfunctions are (f>o — 1/ cosh£ 
and 0i = tanh£. In order to obtain the qualitative behavior of the time dependent solution of 
Eq. ( |152| ), we neglect the contribution from the continuum states with E > 0, only keeping the 
E = —1/2 and E = eigenstates. 

Within this approximation and for sufficiently large times, we obtain 

(159) ^ r ) = ^L + jBtari h£e- T , (r^+oo). 

cosh 4 

where A and B are constant. Returning to original variables, we get 

(160) p(r,t) =p e? (r)|l + Csinh(-^)e 



2R 2 



where p eq is given by Eq. ( f43f ) and C = B/A is a constant. We find that the relaxation time is 
given by t reiax = 2R 2 /a 2 T. 



6 Conclusion 

In this paper, we have studied the Boltzmann-Poisson equation and the Smoluchowski-Poisson 
system in various dimensions of space. Our study shows in particular how the nature of the 
problem changes as we pass from D = 3 to D = 2. We showed that the dimension D = 2 is 
critical in the sense that the results obtained for D > 2 diverge if they are naively extrapolated 
to D — 2. On a physical point of view, the two-dimensional problem differs from the D > 2 
case in two respects: in the 2D case, the central density of the equilibrium state is infinite at 
the critical temperature T c while it is finite at T c in higher dimensions. On the other hand, 
in D = 2, the self-similar collapse results in a Dirac peak which contains a finite fraction of 
mass, while for D > 2, the mass contained in the core tends to zero at the collapse time (but 
a Dirac peak is always formed in the canonical ensemble after t co u as discussed in Appendix 
0). We have also evidenced another characteristic dimension D — 10 at which the nature of 
the problem changes. For D > 10 the classical spiral characterizing isothermal spheres in the 
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physical three-dimensional space disappears. However, since the points in the spiral correspond 
to unstable states, that are therefore unphysical, this transition at D = 10 is not so critical and 
indeed the nature of the self-similar collapse does not show any transition at that dimension. It 
is interesting to note that the dependance of the phase diagram in the (E, T) and (u, v) planes 
with the dimension of space D shows some resemblance with the dependance of the phase 
diagram of confined polytropic spheres with the index n of the polytrope f29fl . An extension of 
our study would be to relax the high friction limit and consider solutions of the Kramers- Poisson 
system and other relaxation equations described in Ref. @. These equations are expected to 
display qualitatively similar behaviors than those described here (i.e. gravitational collapse, 
finite time singularity, self-similar solutions,...) but their study appears to be of considerable 
difficulty since we now need to consider the evolution of the system in phase space. We hope 
to come to that problem in future publications. 
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A Absence of global maximum of free energy in D = 2 

for T <T C = GM/A 

We give a proof similar to the one given in Ref. for D = 3. In two dimensions, we consider 



a homogeneous disk of mass M and radius a at temperature T. It is easy to show that the 
total energy (H) of this disk is 

(161) E = MT + ^GM 2 (\na-^J, 

with the convention $ ~ GMlnr at large distances. According to Eqs. (|6|) and (^), its free 
energy reads 

(162) J = MlnT-Mln^-^ -M-^-flna- 1 
For a — > 0, the free energy behaves like 



2T V 4 



(163) j^2M^l-^lna. 

Therefore, if T < T c = GM/4 the free energy goes to +oo when we contract the system to a 
point. This is sufficient to prove the absence of a global maximum of free energy below T c . This 
also shows the natural tendency (in a thermodynamical sense) of a canonical self-gravitating 
system to collapse to a Dirac peak for T <T C . 



B The case of cold systems (T = 0) 

For T = 0, Eq. ( |38D reduces to 

d64) 

Looking for a self-similar solution of the form (|92|) and imposing the conditions ( |94"D and (|96|), 
we find that the scaling profile satisfies 

(165) xS' + aS= {xS' + DS) S. 

Of course, for T = 0, the exponent a cannot be determined on dimensional grounds, as the 
definition r = a/T/po is n °t relevant anymore. As we will see, a will be determined by 
imposing that the scaling solution is analytic. Eq. ( |165| ) can be readily solved leading to the 
following implicit equation for S: 



(166) - S{x)j = Kx a S(x), 

where K is an integration constant. Now, from the definition of S, we expect a small x expansion 
of the form S(x) = S(0) + ±S"(Q)x 2 + 0(x 4 ), which first implies that 

(167) S(0) = 
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and that (x 2 ) 1 a / D ~ x a , which finally leads to 

2D 



(168) 



n — and K = — ^— ( -\S" 

D + 2 2 V2 1 



D 
D+2 



In terms of the scaling function g(x) associated to the mass profile, Eq. ( |166| ) can be rewritten 



(169) 



2x 



D 



g(x) 



\S»(0)\ 



D + 2 



D + 2 



D+2 
D 



where ^"(O) < is arbitrary. This leads to the exact large x asymptotic behavior 



(170) 



9{x) 



D + 2 VCD + 2)|S"'(0)| 



D 

D+2 D 2 

X~D+2 , 



Moreover, using Eq. (|96|) and (|167|) , we get the exact universal asymptotic behavior of the 
central density 



(171) 



p(0,t) ~ S^itcou-t)- 1 , 



Finally, we note that the implicit equation (|169|) can be written as a parametric set of equations 



(172) 



g{y) 



D + 2 



y, 



x{y) 



y 



D + 2 



S"(0)\y 



D+2 



These results can be obtained by a different, more physical, method. We have indicated in 
Ref. M that, for T = 0, the particles have a deterministic motion with equation 



(173) 



= u = -V$. 

dt 



For a spherically symmetrical system, this can be rewritten 

dr M(r, t) 



(174) 



dt 



where M(r, t) is the mass within r. If a denotes the initial position of the particle located at r 
at time t, we have 

(175) M(r,t) = M(a,0), 
so Eq. (|174|) can be integrated explicitly in 

(176) r D = a D - DM(a,0)t. 
If M(a, 0) behaves like 

(177) M(a,0) = A(a D - Ba D+2 ) + 
close to the origin (which is a regular expansion), then 

(178) M(r, t) = Aa D (l - Ba 2 ), with r D = (1 - DAt)a D + DABa D+2 t. 
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Introducing the collapse time t co u = 1/ DA and considering the limit t — > t co u, we obtain 



D 1 



(179) M(r, t) = — — , with r D = —(t coll - t)a u + Ba 

DtcoU tcoll 



Introducing the scaling variables 



^80) X = T. T^TT' V 



1 



(t coZZ — t)^ 2 ' t co u [(tcoll — t) 1 / 2 _ 

we can put the solution in a parametric form 

(181) M(r,t) = ^(t coll -t) D ^y, with x = (y + Cy^)^ , 
where C is a constant. At the collapse time t = t co u, 

1 Hp" J~) 2D 

(182) M{r,t = t coa ) = —r^+~\ p(r,t = t coll ) = — r~^. 

DC^+i (D + 2)S D C^ 

These results are of course equivalent to those obtained previously. 

We can now use this method to study the evolution of the system for t > t co u (post-collapse 
solution). For t = t co u + St, according to Eqs. (|176|) and (|182[) , the mass contained inside the 
sphere of radius a co u = C~ l l 2 5t^D~ at t = t co u has collapsed at r = 0, resulting in a Dirac peak 
of weight 

(183) M{U ) t) = ^ jm {t-t co ii) D ' 2 . 

Note that in a bounded domain the final collapse to a central Dirac peak containing the whole 
mass occurs in a finite time after t co u. For r > (associated to a > a co u), one has 

1 / 2 

(184) M(r,t) = M(0,t) + r aW- fl ^ 



''coll 



2D 

/ lg5 \ ..!> ... ,.V ( , ( "'-"I! \ P 



Introducing the scaling variables 



r / a , 

(186) x = , y=( ) -1, 

a coll V O-coll , 

we obtain 



|2 
D+2 



1 

2 \ 5 



(187) M(r,t) = M(0,t)(l + y), with x = (1 + y)^ (l - (1 + y)~Ts 

Subtracting the Dirac peak at r = 0, and considering i C 1, for which y ~ ^x D , we find that 
the leading contribution to the mass profile for small r is 

(188) M(r,f)>«^. 
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Hence the density profile does not diverge at r = + for t > t co u. Instead, the density approaches 
the constant value 

(189) p( 0+>i ) = _£_ j 

D4-2 

which decreases with time. The density profile is thus depleted on a scale r ~ a co u ~ 6t m , 
which increases with time. For r ^> a co u, the density profile remains essentially unaffected. 

In principle, the same phenomenon arises for any < T < T c : the density profile obtained at 
tcou ultimately collapses into a central Dirac peak at a time t* > t co u. This solves the apparent 
paradox that the solution at t = t co u has a vanishing central mass and a finite free energy []J. 
In fact, if we allow singular profiles to develop, the evolution continues for t > t co u and the 
Dirac peak with infinite free energy (predicted by statistical mechanics |14|]) is formed during 
the post collapse regime of our Brownian model Q. In practice, degeneracy effects (of quantum 
or dynamical origin) lead to a finite small core of finite density, controlled by the maximum 
allowed degeneracy 



2 As discussed in Sec. 2.1, the results should be different in the microcanonical ensemble. We shall reserve 
for a future communication the full description of the post-collapse regime. 
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